Reconciling Semiclassical and Bohmian Mechanics: 
III. Scattering states for continuous potentials 
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In a previous paper [J. Chem. Phys. 121 4501 (2004)] a unique bipolar decomposition, ^ = 
*i + ^2 was presented for stationary bound states of the one-dimensional Schrodinger equation, 
such that the components $i and ^2 approach their semiclassical WKB analogs in the large action 
limit. The corresponding bipolar quantum trajectories, as denned in the usual Bohmian mechanical 
formulation, are classical-like and well-behaved, even when has many nodes, or is wildly oscillatory. 
A modification for discontinuous potential stationary stattering states was presented in a second 
paper [J. Chem. Phys. !!! !!!! (2005)], whose generalization for continuous potentials is given here. 
The result is an exact quantum scattering methodology using classical trajectories. For additional 
convenience in handling the tunneling case, a constant velocity trajectory version is also developed. 



I. INTRODUCTION 

This paper is the third in a series investigating the use 
of "counter-propagating wave methods" (CPWMs)ii2i2j^ 
for solving the Schrodinger equation exactly. CPWMs 
are a particular variant of the more general multipolar 
decomposition methods, wherein the wavefunction \ff is 
decomposed into two or more components, Thus, for 
a two-term, or "bipolar" decomposition, — ^1 + ^2- 
This rather trivial-seeming procedure can be quite advan- 
tageous in the context of quantum trajectory methods 
(QTMs)^' 5 ' 6 ' 7 ' 8 ' 9 ' 10 ' 11 ' 12 ' 13 i-e. trajectory-based numer- 
ical techniques for performing exact quantum dynamics 
calculations, based on Bohmian mechanics . 14 i 15 i 16 i 17 i 18 i 19 
Conventional QTMs use a single-term or "unipolar" rep- 
resentation of the wavefunction from which all other 
quantities, such as the quantum trajectories themselves, 
are uniquely determined. Multipolar decomposition, on 
the other hand, can lead to radically different QTM be- 
havior for the individual tyj components, owing to the 
fact that the Bohmian equations of motion are non- 
linear 

As applied to wavepacket dynamics for reactive scat- 
tering systems, QTMs suffer from a significant and well- 
known numerical drawback, which to date, precludes a 
completely robust application of these methods. Namely, 
QTMs are numerically unstable in the vicinity of ampli- 
tude nodes and "quasi-nodes" (i.e. rapid oscillations)^ 
owing to singularities in the "quantum potential," Q, 
which together with the classical potential, V, determines 
the quantum trajectories. In the reactive scattering con- 
text, such behavior is always observed, due to interfer- 
ence between the incident and reflected waves. On the 
other hand, if the latter two contributions to the total VP 
were somehow separated, and associated with two differ- 
ent interference-free Vfj components, the node problem 
might well be circumvented. If, in addition, the com- 
ponent field functions were smooth and slowly- varying, 
far fewer QTM trajectories and time-steps might be re- 
quired than for ^/ itself — although the latter could be 
reconstructed at any desired time simply via linear super- 



position of the components. These numerical advantages 
thus provide significant motivation for consideration of 
the bipolar approach. A much more detailed discussion 
may be found in the first two articles of this series, pa- 
per I (Ref. [0) and paper II (Ref. 0). 

The most obvious aspect of any bipolar decomposi- 
tion, including those restricted along the lines of the pre- 



ceding paragraph 
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is that it is not unique. The 



covering function method,—^ for instance, treats W as 
the difference between two very large-amplitude compo- 
nents, thus "diluting" the effects of interference. The 
one-dimensional (ID) CPWM approach^ on the other 
hand, regards the bipolar decomposition 



(1) 



as a superposition of right- and left-traveling counter- 
propagating waves, \&±. For stationary states at least, 
the Eq. ([I]) decomposition is defined such that the *ff± 
components correspond to semiclassical WKB approxi- 
mations, in the large-action limit. In addition to 
providing pedagogical value (semiclassical and Bohmian 
mechanics can not be so reconciled in a unipolar context), 
the semiclassical field functions are typically smooth and 
slowly- varying, i.e. the semiclassical-like CPWM com- 
ponents ^± provide the desirable numerical advantages 
described in the preceding paragraph. 

In paper I (Ref. 0]), a unique CPWM bipolar de- 
composition was determined for bound stationary eigen- 
states of arbitrary ID Hamiltonians. The resultant ^± 
field functions are smooth and interference-free, and ap- 
proach the WKB approximations in the large-action limit 
(within the classically allowed region of space) as desired. 
Moreover, the quantum potentials q± become vanishingly 
small in this limit, so that the bipolar quantum trajecto- 
ries approach classical trajectories, and only a small num- 
ber are required for a numerical propagation, regardless 
of excitation energy, E. In contrast, the unipolar 4" ex- 
hibits motionless trajectories, and an arbitrarily increas- 
ing number of nodes in the large- action/energy limit. Re- 
sults were presented for both the harmonic and Morse 
oscillator potentials. 
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In paper II (Ref. 0), the ID CPWM ideas were mod- 
ified somewhat for stationary scattering states of dis- 
continuous potentials. Although the CPWM decomposi- 
tion of paper I is uniquely specified for any arbitrary ID 
eigenstate — bound or scattering — and in the bound case, 
always satisfies the correspondence principle, the non-L 2 
nature of the scattering states is such that the paper I 
decomposition generally does not satisfy correspondence 
globally. As discussed in Sec. Ill A[ this requires that 
global modifications must be made in order to enable a 
correspondence between and ^>±. These are such 
as to lead to substantial differences between the density 
functions \^> s ±(x)\ 2 and |\I'±(a:)| 2 , although the trajecto- 
ries are identical. Moreover, the resultant \l/± are found 
to correspond to the familiar "incident," "transmitted," 
and "reflected" waves of traditional scattering theory 22 
in the appropriate asymptotic limits. From the time- 
dependent standpoint, reflection was found to be due to 
trajectory hopping from one CPWM component to the 
other, as can be naturally understood using a simple ray 
optics analogy. The method was applied to several ele- 
mentary discontinuous potential systems, including the 
square barrier/ well. 

In the present paper (paper III), the Ref. [2j formu- 
lations are generalized to incorporate both continuous 
and discontinuous potential systems. Once again, an 
analogy is drawn from semiclassical mechanics, albeit a 
"sophisticated" versio n 23 ' 24 i 25 less frequently considered 
(Sec. Ill A[) . As in paper II, a time-dependent method 
based on ray optics is developed for computing station- 
ary states of any desired energy and boundary condi- 
tions. For the most part, the discussion of the preceding 
paragraph still applies, but some additional key points 
should be emphasized. First, the bipolar decomposition 
now provides a sensible definition of "incident," "trans- 
mitted," and "reflected" waves throughout all space, not 
just asymptotically. Second, the explicit hopping of tra- 
jectories from one CPWM component to the other is re- 
placed with a coupling term in the time-evolution equa- 
tions. Third, the trajectories become completely classi- 
cal. Finally, an alternative methodology is also devel- 
oped, based on the use of constant velocity trajectories, 
which can be readily applied to barrier tunneling situ- 
ations. The new methods are found to be remarkably 
efficient, accurate, and robust across a diverse range of 
ID test potentials and system energies (Sec. IIV|) . 

The paper is organized as follows. A derivation and 
discussion of the time-evolution equations for the CPWM 
bipolar components both for classical and constant 
velocity trajectories, are presented in Sec.HTland the Ap- 
pendices. Sec. Mil provides numerical details of the vari- 
ous bipolar algorithms used to compute stationary states. 
Results are presented in Sec. IIVI for four benchmark ap- 
plications: Eckart barrier; square barrier; uphill ramp; 
double- Gaussian barrier. For the first three, these are 
compared with known analytic solutions. Concluding re- 
marks, including prospects for future development, may 
be found in Sec. [Vj 



II. THEORY 

A. Semiclassical Approximations 

Let V(x) be the potential energy for a ID scatter- 
ing system, and E the energy of some stationary state 
^(x) such that E > V(x) for all x. In reality, ^(x) al- 
ways manifests some reflection, even though the energy is 
above the potential barrier. However, basic WKB theory 
predicts zero reflection in this CSS6, clS the classical trajec- 
tories do not turn around. This is also evident from the 
form of the two basic WKB counter-propagating wave 
solutions, 

* s £{x) = r sc {x)e ±is ^ H , (2) 



where 




and s' sc {x) = y/2m[E-V(x)], 



(3) 

m is the mass, F is the invariant flux (paper I), 
and primes denote spatial differentiation. Note that 
since scattering solutions are non-square-integrable, the 
choice of F is arbitrary; however, throughout this pa- 
per we shall adopt the usual left-incident wave normal- 
ization convention, lmx^^oo r sc (x) = 1, so that F = 
lim x _ ) ._ 00 s' sc (x)/m. The positive and negative momen- 
tum functions, p±(x) = ±s' sc (x), specify the classical tra- 
jectories and semiclassical Lagrangian manifolds 1 - ' 26 ' 27 ' 28 
(LMs) associated with the ty±(x) solutions. Clearly, the 
classical trajectories do not change direction, as there are 
no turning points along the real x axis. Therefore, for an 
arbitrary linear combination solution, 

* = a + *^ c + a_* s _ c , (4) 

the entire reflected wave must be due to the $> s °(x) con- 
tribution. The boundary conditions presume that the 
left-traveling wave contribution vanish in the x — > oo 
limit, implying that a_ = and ^S(x) — ty+(x) — i.e., 
there is no reflected wave. 

Various strategies have been developed to deal with 
the above difficulty of the basic WKB method . 23 ' 24 ! 25 
The most common involve analytic continuation into the 
complex plane. Specifically, if no turning point is lo- 
cated along the real axis, then one is found elsewhere 
in the complex plane. The path of the incident wave 
is deformed so as to give rise to a reflected wave upon 
encountering the complex turning point. An analysis of 
the Stokes and anti-Stokes lines 2 ^^ 5 -^^ that emanate 
from the complex turning point and cross the real x axis, 
enables one to effectively recast the above procedure in 
terms of purely real-valued x. The net effect (as inter- 
preted in this paper) is the introduction of coupling from 
= oq. 1 !^ to = a_\l/!f, resulting in an a_ value 
that changes over x, thus yielding meaningful semiclassi- 
cal partial reflection probabilities. 
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Alternatively, there is an approach due to 
Bremme r 23 i 25 i 31 that from the start is formulated 
entirely on the real x axis. This approach is preferred 
for the present purpose, as it provides a common bipolar 
foundation not only for basic WKB and "sophisticated" 
(i.e. capable of predicting partial reflection) semiclassical 
methods, but also for exact quantum scattering applica- 
tions, as will be shown in Sec. Ill Bl Bremmer's idea is to 
model the continuous potential V(x) using a collection 
of discontinuous steps. Solutions are determined for an 
arbitrary step size, which is then made infinitesimally 
small. 

Locally, along any given step, the two "exact" 
Schrodinger solutions (for the model step potential) are 
plane waves, exp(±ipkx), where pk — \J 2m (E — T4 ) , 
and T4 is the (constant) potential value along the fc'th 
step (Appendix A). A viable global solution, ^(x), must 
match boundary conditions appropriately at each of the 
step edges. In reality, a pure positive momentum local 
plane wave in one step would be joined to some linear 
combination of positive and negative momentum plane 
waves in the adjacent step, corresponding to partial re- 
flection and transmission off of the local step. If the re- 
flection contribution is ignored, e.g. so that only positive 
momentum local plane waves are involved, then the re- 
sultant approximate global solution can be easily shown 
to be just ^ B f{x) in the limit of infinitesimal step size. 
Similarly, ^> s ^(x) is obtained from the negative momen- 
tum local plane waves. This approximation thus leads to 
the uncoupled basic WKB solutions of Eqs. and ©. 

The above approach clearly demonstrates how "con- 
tinuous reflection"— — i.e. that arising from continu- 
ous potentials — may be interpreted in more familiar dis- 
continuous terms. It also suggests that reflection re- 
quires coupling between positive and negative momentum 
wavefunction components. This approach has been used 
to develop a sophisticated WKB approximation, simi- 
lar to the Stokes/anti-Stokes approach described above, 
wherein one-way reflection from to , F_ is retained, 
but back reflection is ignored^ On the other hand, if 
all reflection is retained, one can derive a coupled pair 
of first order differential equations for \l/± that exactly 
describes the quantum scattering solution 'F of Eq. |T]). 
Although the ^± components are not themselves sta- 
tionary solutions, they do correspond to the familiar in- 
cident/transmitted/reflected interpretations as discussed 
in paper II. 



Exact Quantum Dynamics Using Classical 
Trajectories 

1. Time- evolution equations 
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FIG. 1: Time-dependent analysis of the ^± coupling, ob- 
tained via trajectory hopping due to local reflection off of 
infinitesimal potential steps. Continuous potential reflection 
is thus described in more pedagogical discontinuous terms. 



using a simple time-dependent approach. The classical 
trajectories and LMs were found to be identical to those 
of the basic WKB approximations, 'Fff (x) — thus satis- 
fying the correspondence principle, and giving rise to 
smooth, well-behaved, interference-free CPWM compo- 
nent functions, $>±. Although the method formally ap- 
plies only to discontinuous potentials, Sec. Ill Al suggests 
that it ought to be extendible to continuous potentials 
as well, simply by modeling the latter using infinitesi- 
mally small steps. This idea is indeed quite straight- 
forward to implement, as discussed in Appendix A and 
Fig. [TJ As in paper II, the result is a trajectory-based 
scattering methodology that is both local and exact in 
its treatment of reflection and transmission, due to the 
time-dependent nature of the approach. In contrast, the 
global character of reflection and transmission in con- 
ventional time-independent exact quantum methods is 
well-recognized. Of the time-independent semiclassical 
approximations, the most sophisticated depend on the 
relative placement of local potential features such as dis- 
continuities and turning pointS ) 23 i 24 ^ although the more 
approximate methods treat such features independently. 

From Appendix A, the hydrodynamic (Lagrangian) 
time-evolution equations for the CPWM bipolar compo- 
nents are found to be 



d^ ± 
~dt 



Zm n 



± 



P_ 

2m 



(5) 



where ±p = ±y / 2m(E — V) are the momenta, which de- 
fine the LM and trajectories. Note that as expected, 
these are identical to those of vE'ff (x), i.e. the trajectories 
are completely classical even though the solutions are ex- 
actly quantum mechanical (Appendix B). All of the other 
desirable properties, as described in the preceding para- 
graph and in paper II, are also found to be true. 



In paper II, an exact quantum, CPWM bipolar de- 
composition scheme for stationary scattering states was 
presented such that a suitable \&± could be constructed 
for a solution 'F with any desired boundary conditions, 



2. Comparison with unipolar Bohmian mechanics 

It is worthwhile to compare Eq. with the 

usual unipolar evolution equations of Bohmian 
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mechanics 4 j 1 6 ' 1 7 i 1 9 Both are exact QTMs, but quantum 
effects are incorporated in very different ways. In 
unipolar Bohmian mechanics, the trajectory evolution 
is determined by the modified potential V + Q, where 
the quantum potential correction Q is responsible for 
all quantum effects. Evaluation of Q requires explicit 
double spatial differentiation of the wavefunction 'J, 
which in turn requires specialized numerical techniques 
fSec. IIIIBl) . Also, Q diverges at nodes, and is otherwise 
numerically unstable in the presence of interference. In 
contrast, the present bipolar QTM scheme avoids inter- 
ference difficulties as desired, owing to the semiclassical 
correspondence. Quantum effects do not manifest in 
the trajectories themselves (as these evolve completely 
classically), and there is no quantum potential in Eq. |5]). 

So where do quantum effects come from in Eq. ([5])? 
As anticipated in Sec. Ill Al these must be due to inter- 
component coupling, i.e. the last term in the equa- 
tion. Note that the amount of coupling is proportional 
to p' oc V'. Thus, the coupling vanishes in the limit that 
V' — > which is reasonable, considering that this is also 
the usual WKB condition and that the basic WKB solu- 
tions are uncoupled. Note that no wavefunction spatial 
derivatives are required in Eq. ([5]) — a decided advantage 
over the unipolar approach. In principle, every trajec- 
tory "splits" into a forward and backward moving tra- 
jectory pair at every point in space and time (paper II). 
The forward-moving trajectory remains on the same LM 
as the source, whereas the backward trajectory "hops" 
onto the other LM. However, unlike the situation in pa- 
per II, this splitting, hopping, and subsequent recombin- 
ing need not be considered explicitly, as it is all implicitly 
dealt with at the differential equation level. Nevertheless, 
conceptually, one may regard trajectory hopping as the 
source of coupling and quantum effects. 

Note that for scattering systems, V — ► in the 
asymptotic limits of x. This implies that there is 
no asymptotic coupling (or trajectory hopping) be- 
tween the two ± components — an essential feature 
from the perspective of computing scattering quanti- 
ties such as global reflection and transmission proba- 
bilities (which according to normalization and bound- 
ary conventions discussed in Sec. Ill Al are obtained re- 
spectively via P rcfi = linx r _ 1 ._ 0o |*_(x)| 2 and P t rans = 
Hindoo \p(x)/p(— x)] \ 1 $> + (x)\ 2 ). Equation ([5|) also en- 
sures that the asymptotic *&±(x, t) solutions are the de- 
sired plane waves, and not some arbitrary linear super- 
position such as a sine wave. Note that in this regard, 
the hydrodynamic time-derivative aspect of Eq. ([5]) is es- 
sential, e.g. the ordinary Schrodinger equation would not 
preclude asymptotic sine wave solutions. 

Although ty± and i&±(x) are identical in the asymp- 
totic limits (apart from the constant scaling factor a±), 
they differ in the interaction region [i.e. a±(x) de- 
pends on x in this region], even though the trajecto- 
ries and LMs are the same throughout. This raises 
some interesting questions vis-a-vis the interpretation of 
standard Bohmian quantities in a coupled bipolar con- 



text, which will be explored more fully in Sec. Ill B 31 
Here, we consider the bipolar quantum potential, which 
in the conventional sense would be obtained via q± = 
— (fi 2 /2m)(r±/r±) with — r± exp(is±/h) and r± and 
s± real. But, it is not clear that such a definition should 
apply in the case where the \f± are coupled and do not in- 
dividually evolve according to the Schrodinger equation. 
Indeed, such a definition would be inconsistent with the 
time evolution of s'_j_, which in any event is itself incon- 
sistent with the classical trajectory evolution (s' ± ^ ±p). 
In this context, it is perhaps more natural to define q± 
such that (V + q±) determines the trajectory evolution. 
According to this definition, q± = for the present bipo- 
lar CPWM formulation, even throughout the interaction 
region. 



3. Additional properties 

As in paper II, the desired stationary state solution, as 
obtained from the Eq. ([5]) evolution equations, is not ob- 
served at all times, but only asymptotically in the large 
t limit. The same ray optics and continuous wave cavity 
ring-down interpretations that apply in paper II also ap- 
ply here. Thus, at t = 0, one starts with a left-incident 
asymptotic plane wave truncated outside the interaction 
region. 

It can be shown (paper II, Appendix B) that as t — > oo, 
within any finite x interval that includes the interaction 
potential, the resultant ^t(x,t) converges exponentially 
quickly to the correct time-dependent stationary state so- 
lution with appropriate x boundary conditions. A proof 
is provided in Appendix B, which, for comparison with 
the time-dependent Schrodinger equation, relies on the 
Eulerian (partial time derivative) version of Eq. |5]). i.e. 



± 



at m 



T^- + Ue-2V) 
Zm n 



± 



P_ 

2m 



+ ■ 



(6) 



Note that Eq. (J6J) involves a single spatial derivative of 
the wavefunction, unlike Eq. ([5]). 

Equation ^ above can be employed to derive a very 
useful flux relationship, 



dp± 
dt 



-j' ± ± ^Re [* + * *_] , 



(7) 



where p± = |\I/±| 2 is the density, and j± — ±(p/m)p± is 
the flux, defined in terms of the actual classical trajectory 
velocities. Taken individually, the dp±/dt do not obey 
continuity (except when p' — 0), implying that the to- 
tal probability for each ty± component is not conserved. 
Together, however, they do satisfy a kind of continuity 
relation, in that <9(p+ + p-)/dt = + j-)', implying 
that the total probability for both and ^_ is con- 
served. This is an important, nontrivial result — quite 
distinct from the usual probability conservation of \& it- 
self. As described in Fig. [21 Eq. ([7|) in effect states that 
dp±/dt is equal to the usual negative flux divergence plus 
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FIG. 2: Schematic indicating the flux, probability, and con- 
tinuity relationships that exist between ^± components for 
the CPWM bipolar decompositions considered in this paper. 
All probability that leaves the component must flow into 
the component. Consequently, J [p+(x) + p~ (x)] dx is 
conserved for all time. 

a coupling term, ±dp cp i/dt = ±(p'/m)Re [*+* rep- 
resenting the rate at which probability flows from one 
CPWM component to the other. 

It should be stated that the classical trajectory CPWM 
bipolar decomposition scheme, in the form described 
above, is not viable for computing stationary eigenstates 
below the potential barrier maximum. Tunneling per se 
is not the issue, since one can in principle apply Eq. ([5]) 
in the tunneling regime using analytic continuation (in a 
manner similar to that applied in paper II) as has been 
confirmed in numerical tests. However, a difficulty arises 
when \p\ — > and \p'\ — > oo, i.e. in the vicinity of the real 
x axis turning points. For discontinuous potentials, this 
only occurs for energies near a piece-wise barrier energy 
Vfc, manifesting as substantially increased propagation 
times (paper II). For continuous potentials, all energies 
below the barrier height exhibit turning points along the 
real x axis, and are therefore problematic. Specifically, 
the p' terms in Eqs. (0 and © lead to numerical insta- 
bilities near the turning points. 



C. Exact Quantum Dynamics Using Constant 
Velocity Trajectories 

1. Motivation 

It should be noted that turning point/caustic issues 
similar to those described in Sec. Ill B 31 are also faced 
by semiclassical methods. Thus, similar remedies may 
presumably be applied here (e.g. complex plane path de- 
formation), although these are not considered further in 
this paper. Instead, we are guided by one of our origi- 
nal motivations, 1 ' 2 to develop methods that avoid such 
semiclassical difficulties altogether. To this end, we turn 
once again to the ray optics interpretation of the bipolar 
approach, as discussed in paper II. 

Note that in classical optics, the "ray" interpretation of 
a given time-evolving electromagnetic field is not neces- 
sarily unique; one has a certain freedom to define rays as 
convenient for a given application^ Consider the case 
of total internal reflection, for instance, at the bound- 



ary between two materials with different indices of re- 
fraction. Quantum mechanically, this corresponds to a 
single-step discontinuous potential at an energy below 
the barrier height (paper II). According to the usual ray 
interpretation, the incident rays refract to the extent that 
they become parallel to the interface, and therefore do 
not penetrate at all into the second medium. This pic- 
ture is physically somewhat incorrect however, in that 
it does not capture the exponentially-damped evanes- 
cent wavei 33 ' 3 Quantum mechanically, this corresponds 
to tunneling into the discontinuous step, which is sim- 
ply not described using the bipolar classical trajectory 
method as presented in Sec. Ill Bl 

On the other hand, a simple, alternative ray interpre- 
tation can be applied to total internal reflection that 
does predict the evanescent wave correctly. According 
to this interpretation, the incident ray velocities remain 
constant across the interface. These rays arc reflected 
within the second medium, giving rise to the evanescent 
wave, and also providing an explanation for the Goos- 
Hanchen phenomenon^ ' 33 ' 35 The quantum-mechanical 
analog would be constant velocity bipolar trajectories, 
for which the incident wave asymptotic velocities remain 
constant throughout the interaction region. Presumably, 
the time evolution equations derived from such a choice of 
trajectories would not depend very sensitively on the en- 
ergy E, even in the vicinity of the barrier height, so that 
tunneling, and classical turning points/caustics would 
pose no special difficulties. 



2. Froman and Froman methodology 

The semiclassical-like CPWM discussed in Sec. Ill Bl 
and the Appendices does not generalize in any straight- 
forward manner for trajectories other than classical. 
However, there is an alternative time-independent for- 
mulation, conceptually similar to the Bremmer (B) ap- 
proach but differing in the details, which does allow for 
such a generalization. This approach is due to Froman 
and Froman (F)^ who use it to define generalized semi- 
classical approximations, although it can also be used 
to derive corresponding exact quantum bipolar decom- 
positions. The guiding principle of the F approach (as 
interpreted here) is that semiclassical solutions should 
satisfy the invariant flux propert y 1 ' 2 ' 24 — i.e., the left side 
of Eq. ([31) with s' sc (x) essentially arbitrary. If s' sc = p is 
chosen classically, then the usual basic WKB solutions 
result for 'i'ff (x). Curiously, however, the corresponding 
exact quantum ^±(x) are not the same as in Sec. Ill Bl 
In the F case, the decomposition of Eqs. |T]) and ^ is 
uniquely defined via the time-independent Schrodinger 
equation and the relation 

= a + (* b ; c )' + a_ (* sc )' (8) 

= -f^ + ?(*+-*-)■ W 
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Thus, in comparing Eq. {8]| to Eq. ([4}, the a± coefficients 
behave in the derivative as if they were x-independent 
constants, though, in fact, they are not (except when 
there is no coupling). 

Although not identical, the F and B bipolar decom- 
positions are somewhat similar in the classical trajec- 
tory case [compare Eq. © to Eq. (|B2|) ] . and both are 
plagued by similar numerical instabilities near turning 
points, as has been confirmed in numerical testing. How- 
ever, the main advantage of the F approach is that it 
can be applied to arbitrary trajectories, p. In partic- 
ular, we choose constant velocity trajectories obtained 
from the left-incident plane waves, i.e. p — \J1mE 
[with lim^^-oo V(x) — 0]. This gives rise to exact quan- 
tum CPWM components \P± that are smooth and well- 
behaved everywhere, even near turning points and bar- 
rier maxima. Note that for constant velocity trajectories, 
Eqs. ([T]) and (|9|) imply that ^± are linear combinations 
of \I> and 



3. Time- evolution equations 



In analogy with Sec. Ill B( the goal is to develop a time- 
dependent method to compute the F constant velocity 
CPWM bipolar decomposition in the large t limit. How- 
ever, since the F construction of the 1 $>± is radically dif- 
ferent from that of B, a substantially different approach 
than that of Appendix A must be used. We have de- 
veloped several different derivations, all of which yield 
the same final results. The simplest strategy is to "work 
backwards" through Appendix B, but using Eq. ([9|) for 
constant velocity trajectories instead of Eq. (|B2[) . How- 
ever, a more pedagogical approach may also be employed, 
as presented below. 

As in Eq. ([8]), the basic idea is to presume that the 
coefficients a± act as constants, but with respect to time 
derivatives rather than spatial derivatives. In particu- 
lar, it seems more natural here to refer to total rather 
than partial time derivatives, given the trajectory-based 
nature of the methodology. Less obvious at this stage 
is the fact that the linear combination — v &-) must 
be used rather than W = fy + + itself, in order to be 
consistent with the time-independent F results. We thus 
obtain the condition 

_ d^ dtf^ 



dt dt dt dt 

iE , p , 

= -—(*+-*_) + —*' 
n m 



(10) 



Together with the time-dependent Schrodinger equation 
(converted to total derivative form), as applied to the 
Eq. ((T|) linear combination, Eq. (JTUJ) above gives rise to a 
unique set of time-evolution equations for d^> ±/dt. In the 
constant velocity trajectory case, for which ty s ±{x, t) oc 

these are found to be 



exp 



i(±y/2mEx - Et)/Ti 



4- Additional properties 

Using a procedure analogous to that described at the 
start of Appendix B, one can show that for station- 
ary states, Eq. (jTTJ) is in fact consistent with the time- 
independent F constant velocity CPWM bipolar decom- 
position. This requires the Eulerian version of Eq. ((11)) , 
i.e. 



at m 



T i(£?-V)* ± -lv* T . (12) 



~dT 



(E - V) ^± 



(ii) 



What is not so clear, however, is whether starting with 
a truncated asymptotic plane wave, one necessarily ap- 
proaches the stationary solution in the large t limit 
(Sec. Ill B 3| paper II). Although this has not yet been 
proven, it is a reasonable assumption, given both the 
counter-propagating trajectory nature of the method and 
flux properties similar to the B classical trajectory case 
(discussed below). Moreover, for all of the numerical ap- 
plications considered in Sec. IIV| exponentially fast con- 
vergence is in fact observed. 

The F constant velocity time-evolution equations 
[Eq. (fTTj) ] offer some decided advantages over the classical 
trajectory approach. Since only energy quantities appear 
on the right hand side, there is no need to resort to an- 
alytic continuation in order to handle tunneling for the 
below-barrier energies. The equations may therefore be 
applied with equal ease throughout the energy spectrum, 
and in fact, the resultant ^± are qualitatively similar 
above, below, and just at the potential barrier maximum 
(Sec. IIV A[) . Another advantage is that the numerical 
propagation does not require spatial differentiation of any 
kind — not even of the potential, V, to determine forces 
driving the trajectories. Furthermore, Eq. (jlip may be 
numerically implemented as is for discontinuous poten- 
tials, just as easily as for continuous potentials, without 
the need for explicit splitting and recombining of trajec- 
tories as in paper II (Sec. HVBj) . 

On the other hand, the F constant velocity approach 
introduces a drawback that the classical trajectory meth- 
ods do not have to contend with when the potential is 
"asympotically asymmetric" — by which we mean simply 
that linij^-oo V(x) = ^ lim^^oo V(x) = Vo, corre- 
sponding e.g. to an exoergic or endoergic chemical re- 
action. Whereas in the x — > — oo limit, the evolution 
equations are uncoupled, this is not true in the x — > oo 
limit if Vo ^ 0, in which case the time-dependent ^f± 
are not expected to converge. Many techniques could be 
used to remedy this situation, e.g. trajectories described 
via smoothly- varying sigmoid (tanh-like) functions rather 
than uniform or classical trajectories. For purposes of 
this paper, we adopt a simpler solution, wherein two 
F constant velocity CPWM bipolar decompositions are 
used, one each for the reactant and product regions. 
Standard boundary condition matching techniques are 
then applied to join these together (Sec. IIII C[) . 

Regarding flux properties, Eq. (|12p can be used to ob- 
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tain 



d P ± ., ^2V T rT 



(13) 



which should be compared with Eq. for the B classical 
trajectory case. Note that the total combined probabil- 
ity for both and \E'_ is conserved here as well, i.e. 
d(p + + p-)/dt = and Fig. [2] still applies. 

The coupling term dp cp \/dt = (2V/h)lm [$ + * is of 
course different from that of Sec. lIIB^l although in both 
cases, the wavefunction inner product cross terms are 
involved. Several additional relations unique to the con- 
stant velocity case may also be easily derived from the 
above equations, e.g. 



dp±/dt = ±dp cp i/dt, 



(14) 



which states that the probability lost by a positive mo- 
mentum trajectory is gained by the negative momentum 
trajectory directly "beneath" it (same x), and 



P+ = P- 



for stationary solutions, 



(15) 



which states that the positive and negative component 
density functions are identical apart from a constant. 
Note that for asymptotically symmetric potentials, and 
the normalization and boundary conditions described in 
Sees. Ill Al and HT~B 21 this constant must be the transmis- 
sion probability itself, Ptrans- 



III. NUMERICAL DETAILS 

In this section, we discuss the numerical details associ- 
ated with implementing the bipolar time-evolution equa- 
tions of Sec. [II] as a practical and robust algorithm for 
computing stationary scattering states. Although any 
boundary conditions may be considered, our emphasis 
shall be on left-incident solutions, for which — > as 
x — > oo . The region of interest is defined via xl < x < 
xr, a region presumed to include the entire potential in- 
teraction to the desired level of accuracy. According to 
the time-dependent ray optics picture developed in pa- 
per II, at the initial time, W = ^ + is a left-incident plane 
wave truncated on the right at x — xl ■ This initial wave 
propagates into the interaction region, wherein it is cou- 
pled to and eventually reaches a stationary state. 
Outside of the interaction region, the 4"± coupling may 
be regarded as effectively zero. Thus, one may compute 
reflection and transmission probabilities, P re fl and Ptrans, 
simply by evaluating ^ + (xn) and ^-(xl) at sufficiently 
large times. 

From a purely numerical point of view, the 
above scheme offers many advantages — e.g., the abil- 
ity to compute specific scattering states without the 
need for complex scalin g 36 ' 37 i 38 or complex absorbing 
potential o 30 ' 32 i 39 ' 40 i 41 ' 42 that would increase the coordi- 
nate range unnecessarily. Moreover, the density of grid 
points may be substantially reduced, owing to the fact 



that the interference- free component functions ty± are 
generally smooth and slowly varying as compared to W 
itself. For the same reason, a larger time-step is also an- 
ticipated. Most importantly however, since q± = 0, there 
is no need to compute on-the-fly numerical spatial deriva- 
tives. With regard to parallel computation, therefore, 
there is no need for trajectories to communicate within a 
LM, although the coupling requires position-dependent 
communication between positive and negative LM tra- 
jectories. 



A. Algorithmic Issues 

In order to solve the time-evolution equations numer- 
ically, two trajectory grids are created, one for each 
bipolar component of the total wavefunction. Hereafter, 
"upper" will be used to describe attributes of the 
component, and "lower" for the component (e.g. 
"lower/upper grid," etc.). On each grid, the correspond- 
ing wavefunction component is spatially discretized at 
the grid point locations and evolved in time. 

For the applications considered here, we have found 
it convenient to modify the ray optics picture somewhat 
from the form described above and in paper II. First, we 
define an initial negative LM grid, even though ^_ itself 
is zero initially. The idea here is that unlike paper II, 
we wish to avoid explicit trajectory hopping. Conse- 
quently, both sets of grid points exist for all time and 
evolve independently of each other (though the compo- 
nent wavefunctions do interact). Second, to avoid un- 
necessary computational effort, no propagation is done 
outside of the interaction region of interest. Instead, at 
uniform time intervals, new ^+ trajectories are fed in at 
x = xl , with an initial value consistent with the pos- 
itive momentum asymptotic plane wave solution. Later, 
as these upper grid trajectories reach x = xr, they are 
discarded. The ^_ situation is similar, except that the 
initial ^^(xr) value is set to zero, and the trajectories 
are discarded as they reach x = xl- 

The most substantial difference from the ray optics pic- 
ture, however, is found in our implementation of the con- 
stant velocity method, for which the grids are distributed 
uniformly throughout the interaction region at all times. 
At the initial time, the upper and lower grids coincide. 
The initial \t + is taken to be the asymptotic plane wave 
solution throughout the interaction region, and the initial 
ty- is set to zero. The above modification, which we call 
the "non-wavefront" approach, is certainly not necessary, 
and is introduced simply for convenience. Calculations 
performed both ways reveal that both converge exponen- 
tially to the correct stationary solution. However, the 
initial convergence of the non-wavefront calculations is 
faster, owing to the fact that the coupling takes effect im- 
mediately. The numerical algorithm is also easier to im- 
plement. Consequently, all constant velocity results pre- 
sented in Sec. IIVI were obtained using the non-wavefront 
code. 
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As per Sec. UH the upper and lower grids move clas- 
sically or with constant velocity, as appropriate, and 
the \l/± contributions evolve in accord with Eq. |5| or 
Eq. Since the grids move in opposite directions, 

they do not remain commensurate over time. Numerical 
interpolation (Sec. lIII B"|) is therefore required to compute 
the coupling contribution from the component wavefunc- 
tion of one grid to the other. As an alternative to the 
above trajectory-based (Lagrangian) approach, we have 
also implemented a fixed-grid (Eulerian) bipolar propa- 
gation algorithm for the constant velocity case. The pri- 
mary advantage of fixed-grid propagation is that the two 
grids remain commensurate for all time, thus avoiding 
the need for interpolation. However, this is achieved at 
the expense of switching to the Eulerian evolution equa- 
tions [Eq. (|12p]. which require explicit numerical spatial 
differentation of the bipolar wavefunction components. 



B. Low Level Issues 

For the numerical propagation of the time-evolution 
equations, two fourth-order explicit integrators were 
considered — the multi-step Adam's/Bashforth, and 
Runge-Kutta methods.— Although, both integrators are 
fourth-order accurate, and require approximately the 
same CPU time, multi-step integrators can not be triv- 
ially implemented for the asymptotically asymmetric po- 
tential case fSec. IIII C|) . Consequently, the basic fourth- 
order Runge-Kutta method was used for all results pre- 
sented herein. In the future, the time integration can 
be made more accurate by using a time-adaptive Runge- 
Kutta approach such as the Cash-Carp method.— For 
most of the calculations performed here, a time-step of 
A = 10 a.u. was found to be sufficient to achieve a com- 
puted accuracy of 10~ 4 or better (Sec. IIV[) . Generally 
speaking however, the time-step issue is quite problem- 
dependent, as it can depend on both grid point velocities 
and spacing. 

As discussed in Sec. IIII A[ the trajectory-based algo- 
rithms require interpolation of both bipolar components 
in order to compute coupling contributions. For all exam- 
ples considered here, these interpolates were calculated 
via a polynomial moving least squares (MLS) method. 4 
In MLS methods, local low-order polynomial fits are cal- 
culated about each grid point via a small stencil of sur- 
rounding neighbor points. To ensure that the fit will 
exactly represent the function values at the grid point 
locations, the stencil size and polynomial order must be 
set equal, thus effectively transforming the MLS method 
into a moving-interpolation method. For all calculations 
performed here, five stencil points were used, correspond- 
ing to a symmetric stencil (for the interior grid points) 
and quartic polynomial interpolation. 

For the fixed-grid algorithms, explicit numerical spatial 
differentiation of both CPWM bipolar components must 
be performed. For the present study, these were calcu- 
lated using centered, fourth-order finite differences^ for 



the interior grid points, and one-sided fourth-order finite 
difference for the boundary grid points at Xl and xr. 



C. Asymptotically Asymmetric Potentials 



As discussed in Sec. Ill C 41 the constant velocity time- 
evolution equations will exhibit coupling in the x — ► oo 
asymptote if V(x) is asymptotically asymmetric, i.e. if 
lirn^oo V{x) = Vq 7^ 0. To remedy this, we con- 
struct two sets of solutions, one for the "reactant" region, 
x < xq, and another for the "product" region, x > xq, 
with xq the dividing point. For the former solutions (de- 
noted via the "L" subscript), Eq. (fTTj) may be used di- 
rectly, with pl — \/2mE. For the product solutions, 
however (denoted "R"), the evolution equations must be 
modified slightly to accommodate the generalized asymp- 
totic potential condition, 



since the positive trajectory momentum is now pa = 
\J2m{E — Vq). It is clear from Eq. (fT6|) that the cou- 
pling vanishes as x — > oo. 

For all propagation times, the wavefunction and its 
first derivative must be continuous across xq. Remark- 
ably, we can treat this boundary condition like an elemen- 
tary plane wave scattering off a step potential (paper II). 
This is achieved via Eq. ([9|), in terms of which the two 
exact conditions become 



(V-Vb)*»F, (16) 



\i/ n 



PL H + = PR [*R+ K-] , (17) 

where = ^l+(xo), etc. In the trajectory-based al- 
gorithm, the left and right incident wavefunction values, 
and are known at any given time, whereas the 

locally transmitted values and ^f L _ are unknown. 

Equation (fT7|) enables one to solve for the two unknowns, 
and thus to continue propagating trajectories through the 
dividing point, Xq. 
The solutions are 



2 PL 



PL + PR 

2PR 
PL + PR 



*L + 



V R- 



PR-PL 



PL - 
PR 



PR 
-PL 



PL +PR 



(18) 



which do indeed correspond to transmission and reflec- 
tion coefficients for stationary states of the discontinu- 
ous step potential (Appendix A and paper II). Numeri- 
cally, the propagation is implemented as follows. Equa- 
tions (fTTj) and (fl"6|) are used until a trajectory reaches Xq, 
at which point it is replaced with a new trajectory on the 
other side of xq via Eq. (|18|) . This requires a "known" 
value from the opposite LM. If the grids are not com- 
mensurate (i.e. trajectories from opposite LMs do not 
arrive at xq at the same time), it is necessary to apply 
extrapolation to the opposite manifold to determine the 
"known" value. 
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FIG. 3: Converged positive component density p+(x) (solid 
line) and negative component density p~(x) (dashed line) for 
the E — 450 cm -1 stationary state of the Vo = 400 cm - 
Eckart barrier system, as obtained using the classical trajec- 
tory CPWM bipolar decomposition. Shaded area indicates 
the Eckart potential V(x). 



IV. RESULTS 

In this section, four different stationary scattering ap- 
plications are considered, all with comparable charac- 
teristic parameters — e.g., the same proton-like mass of 
m = 2000 a.u., barrier height V = 400cm" 1 w 0.0018 
a.u. (note: 10 x smaller than in paper II), and interac- 
tion region, [xl = —3 a.u., xr = 3 a.u.]. All four systems 
were solved using the non-wavefront, F, constant velocity 
trajectory-based method, hereinafter referred to as the 
"constant velocity trajectory" scheme. For the Eckart 
system (Sec. lIV A]) , additional calculations were also per- 
formed using the constant velocity hxed-grid scheme, and 
the wavefront, B, classical trajectory method, now re- 
ferred to as the "classical trajectory" scheme. 

A. The Eckart Barrier 

The first application considered is the symmetric 
Eckart probleni ) 46 i 47 defined via 

V(x) = V sech(ax) 2 1 (19) 

with parameter values Vq = 400 cm" 1 , and a — 3.0 a.u. 
The above potential is plotted on the scale of the inter- 
action region in Fig. [3j 

We have performed numerical calculations for the 
Eckart system using several of the previously described 
algorithms. In the first study, the classical trajectory 
scheme was employed, i.e. the continuous potential ana- 
log of paper II. A time-step of A = 1 a.u. was used, 
and a maximum of 45 trajectories per LM employed at 
any given time. P lcn and Ptrans values were computed 



for various energies E > Vo, to a convergence error of 
10" 4 , and in every case found to match the exact ana- 
lytical results^ to within the same error. A density plot 
of the CPWM bipolar components for a typical solution 
(E = 450 cm" 1 « 0.002 a.u.) is presented in Fig. [3] For 
most energies E, these calculations are roughly as effi- 
cient as the constant velocity calculations described be- 
low. However, the required propagation time does indeed 
increase rapidly as E — > Vq, as expected. In this limit, 
the classical trajectory bipolar decomposition also be- 
comes unstable; the small central peaks evident in Fig. [3] 
become increasingly pronounced, eventually developing 
into singularities. 

The second study is a detailed convergence and ef- 
ficiency comparison between the trajectory and fixed- 
grid implementations of the constant velocity method 
(Sec. IIII Ap . This study also serves as a benchmark for 
establishing the numerical parameter values needed to 
achieve a 10~ 4 accuracy level. Both schemes were ap- 
plied to a calculation of the E = Vq stationary state, i.e. 
the classical "worst-case scenario," for a varying number 
of grid points, N. A time-step of A = 1 a.u. was used for 
the trajectory calculation, and A = 0.1 a.u. in the fixed- 
grid case. Both calculations were propagated to time 
i max = 10, 000 a.u. These parameters were sufficiently 
converged as to ensure that their contribution to the to- 
tal numerical error is insignificant (i.e. 10" 8 or less). In 
Fig. 31 the resultant CPU times on a 2.60 GHz Pentium 
computer are presented. Note that despite the ten-fold 
increase in the number of time-steps for the fixed-grid 
case, the total CPU time is still markedly faster for all 
grid sizes considered. This is due both to the search 
operation (to find the nearest opposite LM trajectories) 
and the interpolation matrix inversions required at each 
time-step by the trajectory scheme. 

In Fig. OJ the log of the computed P ren and Ptrans er- 
rors (relative to known analytical values^) are plotted 
versus N. Across N, the trajectory-based results are 
seen to be much more accurate than the fixed-grid re- 
sults, although the latter are substantially improved by 
increasing the density of grid points. For example, in 
order to achieve 10~ 4 accuracy for both P ren and Ptrans, 
only TV = 25 trajectories were needed, whereas N = 91 
fixed-grid points were required. This is consistent with 
the fixed-grid method requiring explicit spatial differen- 
tiation, introducing an additional source of numerical er- 
ror. Note that for this comparison, the CPU times are 
comparable — 4.3 seconds and 2.7 seconds, respectively. 
The fixed-grid calculation is slightly faster, although this 
situation would be reversed for more accurate calcula- 
tions (e.g. 10" 6 ) and/or higher dimensionalities. The 
performance of both methods is in any event remark- 
able, with the large-iV trajectory case representative of 
the most accurate trajectory-based quantum scattering 
calculations performed to date. Further improvements 
are also possible, however, as discussed in Sec. IIII Bl 

The final study was a constant velocity trajectory 
calculation of P rcn and Ptrans over a wide range of E 
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FIG. 4: CPU time vs. no. of grid points, N, required to prop- 
agate trajectory (empty circles) and fixed-grid (filled circles) 
implementations of the constant velocity CPWM bipolar de- 
composition to a final system time of t ma x = 10, 000 a.u. The 
trajectory and fixed grid time-steps used were A = 1 a.u. and 
A = 0.1 a.u., respectively. 



FIG. 6: Converged positive component density p+(x) (solid 
line) and negative component density p_ (a;) (dashed line) for 
the E — 400 cm -1 stationary state of the Vo = 400 cm -1 
Eckart barrier system, as obtained using the constant veloc- 
ity trajectory CPWM bipolar decomposition. Shaded area 
indicates the Eckart potential V(x). 




FIG. 5: Errors in computed reflection probabilities P rc fl (solid 
lines) and transmission probabilities Ptrans (dashed lines) vs. 
no. of grid points N for the E = 400 cm -1 stationary state of 
the Vo = 400 cm -1 Eckart barrier system, as obtained using 
trajectory and fixed-grid implementations of the constant ve- 
locity CPWM bipolar decomposition. Errors are relative to 
known analytical results (Ref. [47]). 
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FIG. 7: Total wavefunction density p(x) for the E = 400 cm -1 
stationary state of the Vo = 400 cm -1 Eckart barrier system, 
as obtained using the constant velocity trajectory CPWM 
bipolar decomposition (filled circles), and compared with ex- 
act analytical results (solid line). Shaded area indicates the 
Eckart potential Up- 



values — including those above, below, and at the bar- 
rier maximum, E = Vq. The parameters A = 10 a.u., 
tmax — 10,000 a.u., and N = 31 were chosen to cor- 
respond to a target accuracy of 10 -4 . The resulting 
stationary state solution for the E = Vq case is displayed 
in Figs. [5] and [JJ Note the contrast between the smooth, 
slowly- varying CPWM bipolar densities p± in Fig. [6] vs. 
the oscillatory total density p of Fig. [7] Nevertheless, 
as indicated in the latter figure, the numerically recon- 
structed total density agrees with the analytical result^ 



to within the target accuracy of 10 -4 at all grid points. 
Note also that p+ and p_ are identical apart from a con- 
stant, as predicted in Sec. Ill C"4l 

For all 26 energy values considered, the resulting con- 
stant velocity bipolar decompositions are qualitatively 
similar to those presented here for E — Vq. Indeed, in 
no respect do these calculations seem to depend sensi- 
tively on the value of E, unlike the classical trajectory 
case. Consequently, the same numerical parameters as 
described above were used for all energies. The resulting 
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FIG. 8: Computed reflection probabilities P ror ; (filled circles) 
and transmission probabilities Ptrans (empty circles) vs energy 
E for the Vo = 400 cm -1 Eckart barrier system, as obtained 
using the constant velocity trajectory CPWM bipolar decom- 
position, and compared with exact analytical results (solid 
line) . 
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FIG. 9: Converged positive component density p+(x) (solid 
line) and negative component density p-(x) (dashed line) for 
the E — 450 cm -1 stationary state of the Vo = 400 cm - 
square barrier system, as obtained using the constant veloc- 
ity trajectory CPWM bipolar decomposition. Shaded area 
indicates the square barrier potential V(x) 



computed P re fl and Ptrans values are presented in Fig. [8] 
and compared with analytical values. Once again, all 
errors are found to be smaller than 10 -4 . 



B. The Square Barrier 

Although the primary focus of this paper is continu- 
ous potentials, it is worth emphasizing that the methods 
presented here may be applied to discontinuous poten- 
tials as well. As a case in point, we reexamine the square 
barrier system introduced in paper II. Note that the clas- 
sical trajectory scheme for E > Vq would yield results 
identical to those of paper II, as is easily verified from 
Eq. Instead, we focus on the constant velocity tra- 
jectory scheme, which can be applied directly without 
any algorithmic modification. 

For the following square barrier study, we used a 
barrier height of Vo — 400 cm -1 , and barrier edges of 
x\ = —1 and X2 — 1, respectively. As in Sec. IIV A[ 
N = 31 initial grid points were distributed uniformly 
throughout the interaction region, and the time-step was 
defined as A = 10 a.u. The energy was taken to be 
E = 450 cm -1 w 0.002 a.u., i.e. slightly above-barrier. 
The propagation was continued until the computed P ro fl 
and Ptrans values were both converged to less than 10~ 4 , 
which was found to require approximately 1000 time- 
steps. 

Figure [5] is a density plot of the converged CPWM 
bipolar component solutions. As in the Eckart case, the 
results are well-behaved and interference-free through- 
out, although there is a discontinuity in the first spatial 
derivative at the step edges (but not for the total ^ it- 
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FIG. 10: Computed reflection probabilities P rc fl (filled circles) 
and transmission probabilities Ptrans (empty circles) vs energy 
E for the Vo = 400 cm -1 square barrier system, as obtained 
using the constant velocity trajectory CPWM bipolar decom- 
position, and compared with exact analytical results (solid 
line). 



self). For the total density p, the computed values and 
the analytical results are once again in agreement to 10 
or better, at all grid points. 

The calculation described above was repeated for a to- 
tal of 45 different energy values. The resultant converged 
Prefl and Ptrans values are presented in Fig. 1101 as are the 
exact analytical results. Relative to the latter, all com- 
puted probability errors are found to be less than 10~ 4 . 
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FIG. 11: Converged positive component density p+(x) (solid 
line) and negative component density p_ (a;) (dashed line) for 
the E = 500 cm -1 stationary state of the Vb = 400 cm -1 up- 
hill ramp system, as obtained using the constant velocity tra- 
jectory CPWM bipolar decomposition, modified for asymp- 
totically asymmetric potentials. The dividing point is xo = 0. 
Shaded area indicates the uphill ramp potential V(x). 



C. The Uphill Ramp 

The next application considered is an asymptotically 
asymmetric system — the continuous "uphill ramp"— de- 
fined via 



1 



tanh ( — ) 
\2aJ 



(20) 



with the parameters Vq = 400 cm" 1 and a = 0.2 
a.u. This is another scattering system that is ana- 
lytically soluble.— Note that lim^^-oo V{x) — and 
lim^^oo V(x) = Vb, as presumed in Sec. MI Cl 

The propagation was performed using the constant ve- 
locity trajectory scheme described in Sec. IIII C( with the 
dividing point chosen to be io = 0. All other computa- 
tional parameters were identical to those of the previous 
examples, except for the energy, which was chosen to be 
E = 500 cm -1 ps 0.0023 a.u. Again, the propagation was 
continued until both P re fl and Ptrans were converged to 
less than 10 -4 . Figure QT] is a density plot of the re- 
sultant CPWM bipolar component solutions. There is a 
discontinuity at xq = 0, which is to be expected given 
the "step-like" nature of the join (the same behavior is 
observed in paper II). Importantly, however, this discon- 
tinuity does not give rise to any numerical problems, since 
spatial derivatives are not required in the trajectory inte- 
gration (though one must be a bit careful with the inter- 
polation). Away from xo, both density plots are smooth 
and well-behaved. 

Despite the discontinuities in p±, the computed p and 
its first derivative are continuous, owing to the Eq. (|18|) 
relations. This is evident in Fig. [T2J a density plot for 
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FIG. 12: Total wavefunction density p(x) for the E = 
500 cm" 1 stationary state of the Vb = 400 cm -1 uphill ramp 
system, as obtained using the constant velocity trajectory 
CPWM bipolar decomposition (filled circles), and compared 
with exact analytical results (solid line, Ref. [48]). Shaded 
area indicates the uphill ramp potential V(x). 



the converged total $ . As is clear from the plot, no visi- 
ble discrepancies may be observed between the computed 
and analytically exact solutions. 



D. The Double-Gaussian Barrier 

The previous three example systems serve as useful 
benchmarks for testing and evaluating the new bipolar 
methodologies under a wide range of conditions. In par- 
ticular, all three have known analytical solutions, against 
which the computed results may be compared. However, 
we feel it is also worthwhile to consider at least one sys- 
tem which has not previously been solved. One such 
example, which also presents a qualitatively new variety 
of problem, is the symmetric double-Gaussian barrier po- 
tential, 



V(x) = Vb {exp 



-/3(x- x y 



+ exp 



-j3{x + x ) 2 | 



(21) 

For the results presented here, the following parameters 
were used: Vb = 400 cm -1 ; (3 = 9 a.u.; x = 0.75 a.u. 

Once again, the constant velocity trajectory method 
was employed, with the same computational parameters 
as described previously. The energy was chosen to be 
just at the barrier height, E = Vq. Figure [T3l displays 
the resulting converged CPWM bipolar densities. As in 
the previous examples, these are identical apart from a 
constant, and are otherwise smooth and slowly- varying. 
Note that the double-barrier nature of the potential gives 
rise to interesting features in the bipolar densities not 
previously observed. In particular, despite the fact that 
V(x) is changing in the central well region between the 
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FIG. 13: Converged positive component density p+(x) (solid 
line) and negative component density p~(x) (dashed line) 
for the E — 400 cm -1 stationary state of the Vo = 
400 cm -1 double-Gaussian barrier system, as obtained using 
the constant velocity trajectory CPWM bipolar decomposi- 
tion. Shaded area indicates the double-Gaussian potential, 
V(x). 

two barriers, the bipolar densities are nearly flat across 
this region, resulting in a well-defined "reaction interme- 
diate" state between reactants and products. Moreover, 
the density plots enable one to assign quantitative prob- 
ability values for the intermediate state. 

In contrast, the above interpretation and quantitative 
assignments would be very difficult, if not impossible, 
to glean directly from \&(x) itself. This is evident from 
Fig. [141 a density plot of the total \& for the above double- 
Gaussian barrier calculation, obtained via interpolation 
and superposition of the two converged bipolar compo- 
nent solutions, \P±. Note the interference present both 
in the reactant region (due to reflected trajectories) and 
the intermediate region. Thus, not only the intermedi- 
ate probabilities, but also the reflection probability, are 
difficult to read directly from such a plot. 

V. SUMMARY AND CONCLUSIONS 

Scattering applications are of paramount importance 
for chemical reactions, because all reactions may be re- 
garded as scattering events. From a theoretical exact 
quantum perspective, therefore, multichannel scattering 
theory, 22 both time-dependent and time-independent, 
will always play an essential role. At the same time 
however, trajectory-based methods also bring much to 
bear on dynamics, providing great insight into reactive 
processes, vis-a-vis the determination of which trajecto- 
ries make it past the barrier to products vs. those that 
do not. Quantum trajectory methods (QTMs) therefore 
exhibit great potential promise as a chemical dynamics 
tool, combining a trajectory-based description with ex- 




it: (a.u.) 

FIG. 14: Total wavefunction density p(x) (solid line) for the 
E — 400 cm" 1 stationary state of the Vo = 400 cm" 1 double- 
Gaussian barrier system, as obtained using the constant ve- 
locity trajectory CPWM bipolar decomposition. Shaded area 
indicates the double-Gaussian potential, V(x). 



act quantum dynamics. However, all reactive systems 
exhibit interference between the incident and reflected 
(non-reactive) waves, thus causing numerical instability 
problems for conventional unipolar QTMs. CPWM bipo- 
lar decompositions offer a natural means of alleviating 
this interference difficulty, by splitting the reactant re- 
gion ^ into incident \E'+ and reflected "J- components, 
neither of which exhibits interference on its own. More- 
over, this splitting can be extended through the inter- 
action region over to the product region, by which point 
has transformed smoothly into the transmitted wave, 
and \&_ has damped to zero. 

As discussed in Sec. Inland the Appendices, the CPWM 
approach borrows conceptually from semiclassical scat- 
tering methods. Indeed, for the first bipolar decompo- 
sition considered (B, Sec. Ill Bp the bipolar trajectories 
are simply equal to the classical trajectories themselves, 
and the correspondence principle is satisfied in precisely 
the usual WKB limit, V' — > 0. Three features, however, 
contribute to render the present approach fundamentally 
different from basic WKB theory: (1) time-dependent 
formulation; (2) coupling between <3> + and (3) uni- 
versal, local reflection and transmission formulae (see also 
paper II). Point (3) is what determines point (2), i.e. if 
only local transmission were considered without reflec- 
tion, then the coupling would vanish and the basic WKB 
solutions would result. The combination of (1) and (3) 
provides a local physical understanding related to the ray 
optics picture in electromagnetic theory, and also gives 
rise to useful flux relations and numerical algorithms. Re- 
garding the second bipolar decomposition considered, i.e. 
the F, constant velocity scheme, this was motivated by 
practical concerns, but also by an alternative ray optics 
description (Sec. Ill C\i . This can be related to the semi- 
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classical approach of Froman and Froman, and in that 
context, also satisfies a generalized kind of correspon- 
dence principle. 

Note that neither set of evolution equations [Eqs. §5§ 
and (jlip ] involves a quantum potential; instead, all quan- 
tum effects manifest through \I> ± coupling. In both cases, 
the trajectories themselves are not "context-sensitive,"— 
in that they may be computed independently of the ^± 
evolution. Moreover, no spatial differentiation of the 
wavefunction is required, although there may be situa- 
tions where explicit calculation of one spatial derivative 
is numerically advantageous (Sec. IIV A[) . Several other 
numerical modifications have also been introduced for 
convenience fSec. UlT - ^]) . resulting in a shift from an ex- 
act time-dependent Schrodinger interpretation of the dy- 
namics [wherein the exact stationary state is "revealed" 
over time (paper II)] to what may be regarded as more 
of a relaxation approach. Be that as it may, the result- 
ing algorithms offer a remarkably simple, efficient, and 
accurate means of performing reactive scattering calcu- 
lations of all kinds in ID (Sec. lIVp . Indeed, the time-steps 
required for the benchmark molecule- like systems consid- 
ered here are orders of magnitude larger than for typical 
fixed-grid calculations performed at a comparable level of 
accuracy. Moreover, the converged bipolar solution den- 
sity plots render the determination of global reflection 
and transmission probabilities, as well as probabilities 
for reaction intermediate states, quite straightforward. 

In future publications we will continue to generalize 
the methodology described here and in paper II, for the 
type of multidimensional time-dependent wavepacket dy- 
namics relevant to real chemical physics applications. As 
additional motivation for the present work, we now sketch 
how this might be achieved. First, it is necessary to gen- 
eralize the stationary state results of this paper and pa- 
per II for arbitrary time-evolving wavepackets. This is 
done initially for the discontinuous step potential, and 
then generalized for arbitrary continuous potentials in 
a manner similar to Appendix A. In fact, much of the 
groundwork is already laid, in that the time-dependent 
framework has already been introduced. 

The generalization to multidimensional systems is less 
straightforward but can certainly be achieved (such cal- 
culations have already been performed, as will be re- 
ported in a future publication). Conceptually at least, 
many direct chemical reactions can be described using 
a single scattering reaction coordinate, plus additional 
"bound" coordinates. It is natural to consider applying 
the current CPWM bifurcation to the former, and the 
paper I bifurcation to each of the latter. However, the 
total number of wavefunction components would then be 
2 D where D is the number of degrees of freedom. On the 
other hand, for most time-dependent wavepacket calcu- 
lations, node formation in ^ is associated primarily with 
the reaction coordinate itself, due to wavepacket reflec- 
tion off of the reaction profile barrier. Thus, a natural 
approach would be to bifurcate only along the reaction 
coordinate. Only two component wavefunctions result, 



regardless oi D. It remains to be seen whether such a 
procedure will render QTM calculations possible for ac- 
tual molecular systems. Nevertheless, it seems very likely 
that some such bipolar or multipolar approach will go a 
long way towards ameliorating the infamous node prob- 
lem, which has thus far severely limited the effectiveness 
of QTMs in the molecular arena. 
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APPENDIX A: DERIVATION OF BIPOLAR 
TIME-EVOLUTION EQUATIONS 

Following the notation of paper II, Sec. II D 3, we 
presume a discontinuous potential, for which Xk (k = 
1, 2, . . . , I) denote the locations of the I discontinuities. 
The Xk divide configuration space into I + 1 regions, la- 
beled < k < I. In each region fc, the potential has a 
different constant value, Vk- The discontinuous system 
described above may be used to model any continuous 
potential system, V(x), by defining Vk — V{x k n ) [where 
X T = ( x k + Xk+i)/^ is the region midpoint] and taking 
the limit that (xk+i — Xk) — > for all k. 

As the derivation is a time-dependent one, it is conve- 
nient to introduce a small (ultimately differential) time 
increment A, which is then used to determine the Xk 
as follows. Associated with each regi on k is the p osi- 
tive classical momentum value, pk = y2m [E — Vk] [i.e. 
p+(x™)]. The Xk are chosen such that a particle moving 
with momentum pk would traverse the region k in time 
A, i.e. (xt-+i — Xk) = A(pk/m). Consider a trajectory 
which at t = 0, is located at the fc'th region midpoint, x™, 
heading to the right with momentum pk ■ This is clearly 
a positive LM trajectory, carrying a contribution of the 
component wavefunction ^+k = ^+( x = x™ ,t — 0) 
(Fig. [I} . It is presumed that there are also negative LM 
trajectories along which ty- is propagated, but for now 
the emphasis is on 

From t = to t = A/2, the positive LM trajectory 
travels from x = x™ to the discontinuity at x = Xk+i- 
As per paper II, the propagation is that of a plane wave, 
i.e. 



■V+k -* * + (ar fc+ i J A/2) 
fiA 

= * +fe exp - 



2m 
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>\> [ ,rx V [^-[E-2V k ] 



(Al) 



At this point, the trajectory splits into two: one that 
continues in the forward direction along the positive LM, 
transmitting into region fc+1 with momentum Pk+i', the 
other reflected backwards along the negative LM with 
momentum — p k (Fig. [I]). According to paper II Eqs. (17) 
and (18), the trajectory bifurcation introduces a factor of 
2pk/(pk +Pfe+i) hito the the transmitted ^+ wave, and 
{Pk ~Pk+i)/(Pk +Pk+i) into the reflected wave. 

During the remaining time evolution from t — A/2 to 
t = A, the *&_)_ trajectory moves to x™ + 1 [the midpoint of 
the adjacent (k + l)'th region], resulting in an additional 
phase factor analogous to that of Eq. (jAip . The final 
result is 



fiA r 

= exp l—[E-V k - Vfe+i 



(A2) 



Pk + Pk+i 



2p k 



Pk +Pk+i 
Pk - Pk+i 



l + ^E-2V k+1 ] 



,Pk+Pk+i. 

Recall that A = (xk+i — x k )m/p k . In the small A limit, 
Vfc+i -> Vfc, _pfc+i -> _p fc , and -(p fc -p fc+1 )/(x fc+ i - x*) -> 
p'(x) = dp(x)/dx. Substituting into Eq. (|A6|) . and re- 
placing fc subscripts with functions of x = Xk yields: 



(A7) 

2m 



A similar analysis applied to Eqs. (|A4[) and (|A5|) yields: 



eft 



2^ + ^ S - 2F ) 



(A8) 



The reflected trajectory, meanwhile, moves back to the 
original location at x™, resulting in the following for 



* +fc ^*_(x£\A) 



(A3) 



exp [ — [£ - 214] 



Pk - Pk+i 
Pk + Pk+i 



We thus find that *&+ k at time t = contributes 
to both and at time t = A. However, 

these must be combined with similar contributions from 
the initial ^-( k +i) m order to determine the total final 
\I r +(fc+i) and (Fig. [TJ. Following an analysis similar 
to the above, the negative LM contributions are easily 
shown to be the following: 



*_(*+!) -►iMatf+i, A) 
fiA 

= -*_(fc+i) exp I [E - 2V k+1 ] 
*_ (fc+1) ^*_(x^,A) 



(A4) 



= 4< 



/ *A . 

-(fc+i) exp \—\E-V k 



Pk - Pk+i 
Pk + Pk+i 



(A5) 



p k +p k +i 



By adding Eqs. (|A3|) and (|A4[) . we obtain the final ex- 
pression for 1 $ + (x k n +1 , A). Subtracting ^+ k , dividing by 
A, and taking the limit A — > yields the total (hydro- 
dynamic) time derivative, d^> + /dt. We thus obtain 



lim 



lim 

A^O 











VPfc 


+ Pk+l) 



+ -(E-V k -Vk + i) 
h 



APPENDIX B: PROOF THAT BIPOLAR 

STATIONARY SOLUTIONS SATISFY 
TIME-INDEPENDENT SCHRODINGER 
EQUATION 

Consider the stationary solutions of Eq. ©, i.e. 
d^±/dt = —(i/h)E^> ±. Substituting in these expres- 
sions for the time derivatives, and rewriting to obtain 
expressions for the spatial derivatives, yields: 



* '■ = <-£* (B1) 



These equations are identical to those of the bipolar time- 
independent stationary state decomposition described in 
Ref. [25|. Adding the and equations together 
results in 



*' = !<*+ 



(B2) 



Applying spatial differentiation and substituting 
Eq. (|B1|) m to the resulting right hand side yields 
= -(p 2 /h 2 )^. Substitution into the Schrodingcr 
equation then results in 



2m 2m 



W = EV. 



(B3) 



The stationary solutions of Eq. ^ are therefore consis- 
tent with the time-independent Schrodinger equation. 
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